%matplotlib inline
import numpy as np
import matplotlib
import matplotlib.image as image
import matplotlib.pyplot as plt
from skimage.feature import (corner_harris, corner_peaks, plot_matches, BRIEF, match_descriptors)
from skimage.transform import warp, ProjectiveTransform
from skimage.color import rgb2gray
from skimage.measure import ransac
imL = image.imread("images/CMU_left.jpg")
imR = image.imread("images/CMU_right.jpg")
imLgray = rgb2gray(imL)
imRgray = rgb2gray(imR)
# NOTE: corner_peaks and many other feature extraction functions return point coordinates as (y,x), that is (rows,cols)
keypointsL = corner_peaks(corner_harris(imLgray), threshold_rel=0.0005, min_distance=5)
keypointsR = corner_peaks(corner_harris(imRgray), threshold_rel=0.0005, min_distance=5)
extractor = BRIEF()
extractor.extract(imLgray, keypointsL)
keypointsL = keypointsL[extractor.mask]
descriptorsL = extractor.descriptors
extractor.extract(imRgray, keypointsR)
keypointsR = keypointsR[extractor.mask]
descriptorsR = extractor.descriptors
matchesLR = match_descriptors(descriptorsL, descriptorsR, cross_check=True)
print ('the number of matches is {:2d}'.format(matchesLR.shape[0]))
fig = plt.figure(1,figsize = (12, 4))
axA = plt.subplot(111)
plt.gray()
plot_matches(axA, imL, imR, keypointsL, keypointsR, matchesLR) #, matches_color = 'r')
axA.axis('off')
plt.show()
the number of matches is 56
Solution: Let $X$ be the number of matches.
$ 0.95 \geq 1-(1-(1-\frac{N_o}{N})^4)^X \\ X \geq \frac{\ln(0.05)}{\ln(1-(1-\frac{N_o}{N})^)} \\ X \geq \frac{\ln(0.05)}{\ln(1-(1-\frac{32}{53})^4)} = 120 \\ X = 120 $
l_img = keypointsL[matchesLR[:,0]][:,[1, 0]]
r_img = keypointsR[matchesLR[:,1]][:,[1, 0]]
model_robust, inliers = ransac((l_img, r_img),
ProjectiveTransform, min_samples=4, residual_threshold=10, max_trials=20)
matched = matchesLR[inliers]
fig = plt.figure(2,figsize = (24, 8))
axA = plt.subplot(111)
plt.gray()
plot_matches(axA, imL, imR, keypointsL, keypointsR, matched)
axA.axis('off')
plt.show()
fig = plt.figure(3,figsize = (12, 14))
plt.subplot(311)
plt.axis([0, 2*imL.shape[1], imL.shape[0], 0])
plt.imshow(imL)
plt.title("reference frame with the left image")
plt.subplot(312)
imR_warp = warp(imR, model_robust.params, output_shape=(imL.shape[0], 2*imL.shape[1]))
plt.imshow(imR_warp)
plt.title("reference frame with the right image (reprojected)")
im_transform = ProjectiveTransform()
im_transform.estimate(l_img, r_img)
plt.subplot(313)
plt.imshow(warp(imR, im_transform, output_shape=(imL.shape[0], 2*imL.shape[1])))
plt.title("reference frame with the right image (reprojected badly)")
plt.show()
def boundaryDT(image):
r1 = np.arange(image.shape[0])
r2 = np.arange(image.shape[1])
d1 = np.minimum(r1, r1[::-1])
d2 = np.minimum(r2, r2[::-1])
dist = np.minimum.outer(d1, d2)
return np.divide(dist, dist.argmax())
fig = plt.figure(4,figsize = (12, 3))
l_grad = warp(boundaryDT(imL), inverse_map=np.eye(3), output_shape=(imL.shape[0], 2*imL.shape[1]))
plt.subplot(121)
plt.axis([0, 2*imL.shape[1], imL.shape[0], 0])
plt.imshow(l_grad, cmap='gray')
plt.title("Left image DT in Ref. frame (LdtRef)")
imR_grad = warp(boundaryDT(imR), model_robust.params, output_shape=(imL.shape[0], 2*imL.shape[1]))
plt.subplot(122)
plt.imshow(imR_grad, cmap='gray')
plt.title("Right image DT in Ref. frame (RdtRef)")
Text(0.5, 1.0, 'Right image DT in Ref. frame (RdtRef)')
l_grad_adj = np.stack([l_grad, l_grad, l_grad], axis=-1)
r_grad = np.stack([imR_grad, imR_grad, imR_grad], axis=-1)
a_l = np.divide(l_grad_adj, np.add(l_grad_adj, r_grad))
a_r = np.divide(r_grad, np.add(l_grad_adj, r_grad))
l_adj = warp(imL, inverse_map=np.eye(3), output_shape=(imL.shape[0], 2*imL.shape[1]))
l_img_grad = np.multiply(l_adj, a_l)
red_points_img = np.zeros(shape=(imR.shape[0], imR.shape[1], 3))
for point in keypointsR[matched[:,1]]:
red_points_img[point[0]:(point[0]+5), point[1]:point[1]+5] = [255, 0, 0]
red_points_warp = warp(red_points_img, inverse_map=model_robust, output_shape=(imL.shape[0], imL.shape[1]*2))
l_points = keypointsL[matched[:,0]][:,[1, 0]]
fig = plt.figure(5,figsize = (12, 14))
plt.subplot(311)
plt.axis([0, 2*imL.shape[1], imL.shape[0], 0])
plt.imshow(l_img_grad)
plt.title("alpha based on DT (for RransacRef)")
r_img_grad = np.multiply(imR_warp, a_r)
plt.subplot(312)
plt.imshow(r_img_grad)
plt.title("alpha based on DT (for RransacRef)")
pan = np.add(np.add(l_img_grad, r_img_grad), red_points_warp)
axD = plt.subplot(313)
plt.imshow(pan)
plt.plot(l_points[:, 0], l_points[:, 1], '.b')
plt.title("Panorama")
plt.show()
C:\Users\HP\AppData\Local\Temp\ipykernel_20572\186063864.py:4: RuntimeWarning: invalid value encountered in divide a_l = np.divide(l_grad_adj, np.add(l_grad_adj, r_grad)) C:\Users\HP\AppData\Local\Temp\ipykernel_20572\186063864.py:5: RuntimeWarning: invalid value encountered in divide a_r = np.divide(r_grad, np.add(l_grad_adj, r_grad)) Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).